Embedding dissipation and decoherence in unitary evolution schemes 
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Dissipation and decoherence, and the evolution from pure to mixed states in quantum physics are 
handled through master equations for the density matrix. By embedding elements of this matrix in 
a higher-dimensional Liouville-Bloch equation, the methods of unitary integration are adapted to 
solve for the density matrix as a function of time, including the non-unitary effects of dissipation and 
decoherence. The input requires only solutions of classical, initial value time-dependent equations. 
Results are illustrated for a damped, driven two-level system. 
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The study of open quantum systems is of widespread 
interest across different areas of physics particularly in 
the irreversible processes of dissipation and decoherence 
afforded by coupling to an external reservoir or environ- 
ment. Quantum optics is replete with such studies for 
optical bistability, resonance fluorescence, and the gen- 
eral evolution from pure to mixed states, often consid- 
ered through damped, driven two-level atoms (l). Cou- 
pled quantum wells in a wider context and the study of 
quantum Brownian motion, dissipation and fluctuations 
have also received much attention [Q. Application of 
such considerations to "quantum non-demolition" in the 
emerging field of laser-interferometric gravitational wave 
detection, and of quantum noise and decoherence in the 
field of quantum computation, add to the importance of 
this subject. Finally, this evolution from pure to mixed 
states is at the heart of the problem of measurement in 
quantum theory ||. 

On the other hand, unitary integration schemes for 
the evolution operator of time-dependent Hamiltoni- 
ans, when available, are powerful because they pre- 
serve invariants and are stable, also in numerical ap- 
plication. In this Letter, we present a general proce- 
dure and illustrate with an example how to preserve 
most of these advantages even while working with sys- 
tems exhibiting dissipation and decoherence. There 
are two key steps. First, the n-dimensional Liouville- 
von Neumann-Lindblad (LvNL) equation containing dis- 
sipation and decoherence is embedded in a (n 2 — 1)- 
dimensional Liouville-Bloch form with a non-Hermitian 
Hamiltonian. Second, this Liouville-Bloch equation is 
handled by a "unitary integration" procedure that has 
been described in recent years Q ^ |(| wherein the evo- 
lution operator is written as a product of exponentials, 
each exponent involving an element of a closed Lie alge- 
bra of operators together with a multiplicative classical 
function of time. With all the non-commutativity han- 
dled analytically, the entire problem is reduced to solving 
coupled, first-order differential equations for this set of 
classical functions. In many cases, this set reduces to a 
single non-trivial Riccati (first order, quadratically non- 
linear) equation for one of the classical functions, all the 



rest then obtained through trivial quadratures || . All of 
the above features remain valid even when the Hamilto- 
nian is non-Hermitian and the evolution non-unitary. 

Two other papers share our aims in setting the pas- 
sage from pure to mixed states in a unitary evolution 
scheme but they proceed differently. One deals with weak 
dissipation, handling the Hermitian part of the LvNL 
equation through unitary integration and the dissipative 
terms through conventional integrators @. Because of 
their focus on numerical integration, both these han- 
dlings are for small time steps whereas we aim for in- 
tegration over arbitrary, finite t. Another work Q in- 
troduces a novel "square root operator" of the density 
matrix and an associated n 2 -dimensional Hilbert space, 
along with additional constraints that are not in conven- 
tional quantum mechanics. Our embedding in a higher 
dimensional space does not introduce any new elements 
beyond those already in the density matrix. After sub- 
mitting our Letter, we have learnt of another work that 
solves master equations by invoking an "auxiliary" n 2 - 
dimensional Hilbert space 

We begin with the master equation for the density 
matrix p, sometimes called the Liouville-von Neumann- 
Lindblad equation m [| ||, 

k 

= [H, P ] - \(£ (L\L kP + P L{L k - 2L kP L\) ,(1) 

k 

where an over-dot denotes differentiation with respect to 
time and h has been set equal to unity, H is a Hermi- 
tian Hamiltonian, and the second term on the right-hand 
side is the "Liouvillian super-operator" describing cou- 
pling to the environment and the resulting irreversibil- 
ities of dissipation and decoherence. The above form 
in the Markov approximation with an explicitly trace- 
less right-hand side guarantees conservation of Tr(p) and 
positivity of the probabilities. For a more mathematical 
description in terms of so-called "dynamical semigroups," 
we refer to @ |TTJ . 
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Our aim in this paper is to solve Eq. |l]) for fairly gen- 
eral time-dependences of H and the L's contained in it, 
while keeping as closely as possible to the unitary inte- 
gration that applies in the absence of the super-operator. 
This method [|, |(| has been developed when H(t) is a 
sum of terms, each of which involves a time-independent 
operator multiplying a classical function of time. In such 
a case, without any recourse to time-ordered Dyson ex- 
pansions, one can solve for the evolution operator U(t) 
satisfying 



where one convenient choice for the (n 2 



iU(t) = H(t)U(t), U(0)=1, 
by writing U(t) as a product 

U(t) = '[[exp[-in J (t)A j ], 



(2) 



(3) 



where Aj are the operators contained in Hit) together 
with a sequence of other operators formed out of their 
mutual commutators in a successive fashion. If this set 
forms a closed algebra under commutation, then upon 
substitution, Eq. (|^) can be shown to satisfy Eq. (^) 
through repeated application of the Baker-Campbell- 
Hausdorff (B-C-H) identity [4,6]. This results in a well 
defined set of coupled first-order, generally nonlinear, 
equations for the functions Pj(t). Thereby the quantal 
problem is reduced to the classical one of solving this set 
of equations, following which p(t) is obtained as 



p(t) = U{t)p{Q)U\t). 



(4) 



In extending this procedure to non-unitary evolution, 
if we were to retain only the first two terms in the 
superoperator, it is simple to extend Eq. Q by us- 
ing two different products Ux,(t) and Ur(£) so that 
p(t) = Uttf) p{0) U R (t), with correspondingly different 
functions PLj(t) and PRj(t) in Eq. (||). Once again, upon 
calculating ip with such a form, the B-C-H identity can 
be used to get a well-defined set of equations for the pl 
and p,R. However, the last term in the superoperator in 
Eq. ([!]), wherein p(t) occurs between operators multi- 
plying it both on the right and from the left, no longer 
permits easy generalization. Note that this last term is 
the so-called "quantum jump" in interpretations of the 
LvNL equation as conventional continuous evolutioni, al- 
beit with a non-Hermitian Hamiltonian, plus a jump [fl2"[ . 

For the full master equation, we proceed by separating 
the invariant Tr(p) from the n 2 elements pij(t). Eq. ([j]) 
then reduces for the remaining n 2 — 1 elements to the 
Liouvillc-Bloch form 



r\ is pn 



2,3,. 



Pij + Pji, 



— 1) elements of 
Pij - Pji, i > j- 



The first (n — 1) of these describe the diagonal elements 
of the density matrix, the other (n 2 — n) i ^ j, describe, 
respectively, in-phase dispersive and out-of-phase absorp- 
tive components of polarization. Even though C may not 
be Hermitian, the form of Eq. (^|) is now the same as in 
Eq. (Q) with all operators to the left of r\ so that the 
same procedure of a product exponential form for r](t) 
as in Eq. (H) can be carried out now in the (n 2 - 1)- 
dimensional space. Thereby, the LvNL equation for p has 
been embedded in a higher-dimensional Liouville-Bloch 
equation. While invariants are no longer preserved with 
C non-Hermitian, the advantages of exponential factors, 
with all operator aspects handled analytically and only 
classical time-dependent equations to solve, still remain. 

One immediate consequence is worth noting. If the 
operators Lk in Eq. (Q) are such that C in Eq. (0) in- 
volves imaginary elements and, consequently, r\ decays 
asymptotically, r\(t — > oo) — *■ 0, then all coherences van- 
ish (off-diagonal py) and all diagonal pa become equal, 
pu(t -»■ oo) (l/n)Tr(p(0)). Tr(p 2 ) on the other hand, 
decreases asymptotically to (1/n) of its initial value. A 
specific n — 2 illustration will be given below of this 
rather general conclusion. 

To demonstrate this method, we turn now to a series 
of recent papers |l3| that discussed phase coherences and 
transitions in a periodically driven two-level system with 
a single L in Eq. (m: 



H = -e(t)a z + J(t x , L = Vfcr z , pij(0) = SySn. (6) 

Applying our procedure, we have pn (t) + P22 (*) = 1> an d 
Eq. (ph for the three remaining elements takes the form 
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if,{t) = £(t)T){t), 



(5) 



To solve this as a product of exponentials, we need the 
eight operators of an SU(3) algebra. Instead, we illus- 
trate first a simplified variant of Eq. (|^) as our model, 
with a symmetric choice for the Lk involving all three 
Pauli matrices, that is, Lk = ^T/2ak- This modifies 
Eq. (Q) to introduce also a (—iT) in the third diagonal 
element of the matrix. With the matrix then expressible 
as 

C = -iTl - e{t)A z + 2JA X , (8) 

where A x , A y , A z are the operators of angular momen- 
tum in a representation 
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(9) 



the closed Lie algebra of these three suffices to solve 
Eq. (Q) by our unitary integration procedure. Since this 
procedure rests only on the commutators between Aj , we 
can use any representation of them as is convenient. We 
exploit this in choosing Eq. ([)]) so that C involves the Aj 
only linearly. Although, for comparison with |fl3|| , only 
e in Eq. (|^) is a function of time, we note that every- 
thing that follows applies also to more general time de- 
pendences of J and r and inclusion of a time-dependent 
term in A y as well. We also note that reduction of the 
term involving T in Eq. (^) to a unit operator reflects a 
general sum rule in any dimension n. When k in Eq. (|l|) 
runs over all n 2 linearly independent operators, that sum- 
mation reduces to 2(npij — Sij). 

The first term in Eq. (JsJ) leads to a trivial factor 
exp(— Tt) and the remaining Hcrmitian part of C has 
been solved before H: 



rj(t) = cxp[— Ft] exp[— ip + (t)A + ] 

x exp[-i^_(t)yl_] cxp[~ifi(t)A z ]^(0), (10) 



with A-. 



A x ±iA y , 77(0) = (0,0,1), and 
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FIG. 1: p22(t) for an oscillating driving field with J/lu — 3, 
A/uj — 45, and damping values (a) F/lo — 0, (b) F/cu = 0.35, 
and (c) F/lu — 5. 



[i + - ie(t)p + - J(l+fi 2 + ) 
(1 = 2iJfi + - e(t), 
/i_ — = J, 



= 0, 



Mo) 



(lla) 
(lib) 
0. (11c) 



The first of these equations, involving p+(t) alone in Ric- 
cati form, is the only non-trivial member of this set. So- 
lutions give through Eq. (|l0|), 



- ^ + ^exp(-rt)[l 

P22 {t) - i[i-exp(-rt)] + 

p 12 (t) = ip,-(t)exp(-Ft), 
pn{t) = ip+(t)[p,+ (t)p-(t) 



-2/i+(t)yi_(t)], 
A*+(*)M-(t)exp(-rt), 



1] exp(-H). 



(12) 



These are general solutions, valid for any time. The co- 
herences vanish asymptotically and p\\ and p 22 attain 
the value \ as t — > 00. While Tr(p) remains always at 
unity, Tr(p 2 ) decreases to (1/2). The above assumed as 
initial state the pure state with pn(0) = 1 the only non- 
zero elememt, but a wider choice also leads to the same 



final result. Simple numerical integration of Eq. (lla 



for an oscillating driving field e(t) — Acos(uit) are shown 
in Figs. 1 and 2 for various values of the parameters 
(uj, J, A,F) . They are in agreement with In Fig. 

2(c), we also record the time evolution of the entropy, 
S = —Tr(p In p). The value of F governs the rate of rise 
as S increases monotonically from to its asymptotic 
limit of In 2. 

We already noted from the 3x3 matrix structure of 
Eq. (Q) that for the most general H and L in Eq. (^), 
a product of eight exponential operators always provides 
the requisite 77(f). As another illustration of a smaller 
set sufficing, when only three of the four linearly inde- 
pendent matrices are included in Lk, an additional in- 
homogeneous term in the column vector — tT(0, 0, 1) ap- 
pears on the right-hand side of Eq. (^), with L again as in 
Eq. (||). For J = 0, this is easily solved, diagonal and off- 
diagonal elements decoupling, and gives the result that a 
mixed state evolves to the pure state (1,0). 

The reduction in the number of exponential factors 
required is a generic feature, whenever C in Eq. (|5|) 
involves only the elements of a sub-algebra of the full 
algebra of SU(n 2 — 1). Thus, in the n — 2 exam- 
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pies considered above, the existence of SU(2) subalge- 
bras allows solutions with just three exponential op- 
erators in Eq. (|l0|). Denoting the eight operators of 
SU(3) by Oi,i — 1 — 8, with one choice for them being 
(A z , A+ , A_ , A\ , A\ , A 2 _ , A+ A 3 + A 3 A + , A_ A 3 + A 3 A- ) , 
there are several triplets that close under commutation. 
These include the familiar i = (1,2,3) as in Eq. (||) but 
also many others such as (1,5,6) and (1,7,8). There are 
also sub-algebras involving four (for example, (1,2,5,7) 
and (1,3,6,8)) and five elements (examples: (1,2,4,5,7) 
and (1,3,4,6,8)) in which case four or five exponential 
factors, respectively, would suffice for our solution in 
Eq. (ficil). As n increases, although the total number of 
operators n 2 — 1 grows rapidly, once again, C may involve 
only the operators of sub-algebras, SU(n 2 — 1) contain- 
ing many sub-algebras of lower order all the way down to 
SU(2) with just three operators. Indeed, with increasing 
n, there are many more such sub-algebras so that very 
often the H and Lk may afford reduction of the number 
of exponentials in our procedure to a small number. 

In summary, an n-dimcnsional LvNL equation describ- 
ing dissipation and dccohcrcncc (or, alternatively, con- 
tinuous evolution plus a quantum jump) of the density 
matrix p(t) is first embedded into an (n 2 — l)-dimensional 
Liouville-Bloch equation for diagonal and off-diagonal 
combinations r)(t) of pit). A unitary integration scheme 
is then applied to this form of the equation, with rjit) 
expressed as a product of exponentials involving a lim- 
ited, finite number of factors and operators, often just the 
three of angular momentum. Through this procedure, all 
elements of pit) are obtained in terms of solution of a sin- 
gle Riccati equation for a classical function together with 
ordinary multiplication and integration. 

We thank Drs. Dana Browne and Lai Him Chan for 
suggesting we follow the entropy of evolution. 
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